


CHARACTERIZING TWO-TIMESCALE NONLINEAR DYNAMICS 
USING FINITE-TIME LYAPUNOV EXPONENTS AND VECTORS* 

K. D. MEASEt U. TOPCUt AND E. AYKUTLUGt 

Abstract. Finite-time Lyapunov exponents and vectors are used to define and diagnose boundary- 
layer type, two-timescale behavior and to construct a method for determining an associated slow 
manifold when one exists. Two-timescale behavior is characterized by a slow- fast splitting of the 
tangent bundle for the state-space. The slow-fast splitting defined by finite-time Lyapunov expo- 
nents and vectors is interpreted in relation to the asymptotic theory of partially hyperbolic sets. The 
method of determining a slow manifold developed in this paper is potentially more accurate than 
op an existing approach that is based on local eigenvalues and eigenvectors, at the expense of more 

^^ computation, and is more generally applicable than approaches, such as the singular perturbation 

Q_^ method, that require a special coordinate representation. The approach is illustrated via several 

Cn examples. 

^ Key words. Nonlinear dynamics; Multiple timescales; Slow manifold; Lyapunov exponents and 

^~5 vectors 

^^ 1. Introduction. When a finite-dimensional nonlinear time-invariant (NTI) dy- 

^O namical system evolves on multiple timescales, reduced-order analysis may be possible. 

We consider only two timescales in this paper, referred to as fast and slow, but the 

discussion and results are relevant to systems with more than two timescales as well. 

Multiple timescales may induce geometric structure in the flow on the state-space. If 

Cd the system differential equations can be expressed in terms of coordinates adapted to 

H this structure, the system can be decomposed into lower-order subsystems to simplify 

'— ^ the analysis of the dynamical behavior, analogous to modal decomposition for linear 

^^ time-invariant (LTI) dynamical systems. 

^ As a motivating example, consider the optimal (minimizing a weighted sum of 

^\ time and fuel consumption) flight of an aircraft between distant locations. The first- 

C^ order necessary conditions for the optimal solution take the form of a boundary- value 

^^ problem for a Hamiltonian system. The solution may have a "take-off/ cruise/landing 

^^ structure" [24 . The aircraft spends most of the time in the cruise phase flying on 

r ■ a hyperbolic slow manifold in the state-space for the Hamiltonian system where the 

^^ most time/fuel eflicient flight is possible. On the slow manifold, the aircraft travels 

/— s from the vicinity of the (longitude, latitude) of the take-off location to the vicinity of 

the (longitude, latitude) of the landing location. However, this cruise segment does 
not satisfy all the boundary conditions; for example, the altitude may be 35,000 feet. 
The take-off and landing phases are transitions to and from the cruise segment on 
^ the slow manifold, involving some fast behavior in comparison to the behavior on the 

slow manifold. For on-board guidance purposes, it would signiflcantly simplify the 
problem to treat the cruise guidance (guidance on the slow manifold) and the take-off 
and landing guidance (guidance in the fast boundary-layers) separately, reducing the 
order of the relevant dynamics in each case and reducing the numerical sensitivity 
in determining the optimal flight path. This conceptual example has been used be- 
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cause it is easy to describe. For an actual, detailed example see [41] in which such a 
reduced-order approach is proposed for guiding an air-breathing launch vehicle in a 
near- minimum-fuel manner. Minimum- fuel path planning is carried out on the slow 
manifold and, because of the order reduction and reduced sensitivity, this is expected 
to be feasible for near-real-time, onboard computation. A suboptimal feedback con- 
trol law is then used to track the slow manifold trajectory. Because most of the flight 
occurs on the slow manifold, the performance is near-optimal. In |41j . the order re- 
duction and slow manifold characterization were achieved via the analytical singular 
perturbation method, with the geometric singular perturbation theory f51, described 
later, providing insight. 

The analytical singular perturbation method [24l [34] starts with the "standard 
form" of a singularly perturbed dynamical system 



y = ep(y,z,e), 
z = q(y,z,e), 



(1.1) 



where the system state is (y,z), e is a small parameter, and both p and q are of 
0(1) with respect to the e — > limit. The 'dot' over a variable denotes differentiation 
with respect to time t. With the standard form, the terms slow and fast are due to 
the properties y — 0{e) and z = 0(1). For e = 0, there is a manifold of equilibria 
given by £ = {(y^z) £ K" : q(y,z,0) = 0}. Assuming that the coordinates have 
been chosen so that 9q/9z is nonsingular, there exists locally a function ■0(y, e), such 
that q(y, ip{y, 0), 0) — 0, with which the slow manifold can be represented as a graph 
S{e) = {(y,z) G M" : z = ip{y,e)}. The function tp satisfies the partial differential 
equation (PDE) 

(l{y,2.,e) =e--{y,e)p{y,z,e), (1.2) 

oy 

and can be computed as an asymptotic expansion in s. Moreover there is a systematic 
procedure for constructing a solution for particular boundary conditions by matching 
asymptotic expansions for the segment on the slow manifold and fast boundary layer 
segments, though determining terms beyond zeroth-order may be difficult, yet needed 
for sufficient accuracy. A variation on this approach is to pose the PDEs for the slow 
manifold directly. This approach has been developed and applied in the chemical 
kinetics context |5l|36]. However, the model required for the successful application of 
this approach must either be in singularly perturbed form or in coordinates for which 
it is known how to separate them into dependent and independent variables for the 
representation of the slow manifold as a graph, which implies knowing the dimension 
of the slow manifold as well as its orientation. Also the solution of the PDEs is not 
always feasible [S]. 

An important step toward greater generality in modeling and characterizing two- 
timescale behavior was taken by Fenichel [8j|20j. He developed a geometric (coordinate- 
free) singular perturbation theory for two-timescale behavior starting from a given 
one-parameter family of vector fields x — fe(x), e a small constant parameter, avoid- 



ing the a priori need for the standard form representation (1.1 1. He assumed that 
£ = {x. E M" : fo(x) = 0} is a normally hyperbolic manifold of equilibria (fixed points) 
and showed that this manifold persists in perturbed form for small nonzero values of 
e. For nonzero e, the perturbed manifold, denoted by S{e), such that 5(0) = £, is 
no longer composed of equilibria; however it is invariant under fe(x) and for small s 
the motion on S{s) is slower than the transverse motion, so it is referred to as a slow 
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invariant manifold. Conveniently eigen-analysis (meaning the use of matrix eigen- 
values and eigenvectors) of Z?fo(x) = ^(x) provides the timescales and state-space 
structure near 5(0) = 5 in the extended space M" x M, where the extra dimension 
is for £. Associated with each x g 5(0) are center, stable and unstable subspaces of 
the (n + l)-dimensional tangent space. The center manifold theorem [15] connects 
the center, stable, and unstable subspaces to the corresponding nonlinear invariant 
manifolds. For a particular value of e, the slow manifold S{e) is a slice of the center 
manifold. The left side of Fig. (1.1) depicts the geometry for n = 2 in the extended 
space. 




■•£22 = ?(«/, 2, £2) 



r = Ux) 




Fig. 1.1. Fenichel's geometric perspective of a two-timescale system. 



Starting with a general system x = f (x), as depicted for n — 2 on the right side of 
Fig. (1.1), one strategy is thus to suspect two-timescale behavior with a slow invari- 



ant manifold S and seek a coordinate transformation x — h(y,z,e), where y € 



and z e 



ith 



J = 



for which the dynamics take the standard form of 



(1.1 1 with the requisite properties for the application of the singular perturbation 



method. The singular perturbation method has indeed found utility in a number of 
fields [24l [33l [34j • However, its applicability is limited, for lack of a systematic method 
of diagnosing two-timescale behavior and generating the required standard form ( 1 . 1 1 



representation, whose special coordinates are compatible with representing the slow 
manifold as a graph. The air-breathing launch vehicle application mentioned above re- 
quired hypothesizing the two-timescale behavior, choosing an appropriate coordinate 
representation, and artificially introducing a small parameter into the representation. 
The question arises as to whether two-timescale behavior can be defined inde- 
pendently from a singularly perturbed dynamic model, and if so, whether there is a 
methodology by which it can be diagnosed and exploited for reduced-order modeling 
and analysis. Timescales, state-space structure, and reduced-order analysis are ap- 
plicable to nonlinear dynamical systems of the general form x = f (x) on M", but are 
significantly more challenging to define and implement. One approach is to analyze 
the linear variational equation v = I?f (x)v along orbits of the nonlinear system to 
identify the timescales and the associated tangent space structure. The tangent space 
structure is then "transferred" to the local manifold structure of the nonlinear fiow 
in the state-space. When considering asymptotic behavior, the usual tangent space 
splitting m 122] is into stable, center, and unstable subspaces; these subspaces are 
transferred to the corresponding manifolds. For finite-time behavior, we want to dis- 
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tinguish behavior that either dies out quickly in forward or backward time, relative 
to the characteristic time interval of interest T^, from slower behavior. If the "orbit" 
under investigation is an equilibrium point Xg, the variational system is LTI, so eigen- 
analysis of Df (xe) is applicable and can yield an invariant splitting of the tangent 
space at Xg, denoted by Tx^M" — E^'^^Xe) i?*(xe) © E^'^iTCe), into fast contracting, 
slow, and fast expanding linear subspaces. These subspaces can be transferred to fast 
contracting, slow, and fast expanding manifolds for the nonlinear flow near Xg |15| . 
For a periodic orbit, eigen-analysis can be applied to the corresponding linearized 
period-one map to determine the tangent space structure, and then this structure 
transferred to the manifold structure of the nonlinear flow in the neighborhood of 
the periodic orbit. Away from equilibria and periodic orbits, the characterization of 
timescales and associated state-space structure is more difficult. 

An existing approach to analyzing two (or more) timescale behavior for a gen- 
eral coordinate representation x — f (x) away from equilibria and periodic orbits - 
that does not rely on characterizing it as a singular perturbation - is the intrinsic 
low-dimensional manifold (ILDM) method [33], developed in the context of chemi- 
cal kinetics. Via eigen-analysis of the system matrix Di{x.) for the linear variational 
dynamics at points x in the state-space region of interest, the tangent space is split 
into slow and fast subspaces: TxM" = _B*(x) © E^{x), where the "hat" denotes that 
these subspaces are approximations, as clarifled later, of the slow and fast subspaces; 
in particular these subspaces are not invariant with respect to the linear variational 
dynamics. If a slow invariant manifold S exists, the invariance of the slow manifold 
with respect to the flow of x = f (x), and the fact that the flow on the slow manifold is 
slow, dictate that at points on S, f (x) should lie in the slow subspace, or equivalently 
that f (x) should be orthogonal to all the vectors in the orthogonal complement to 
the slow subspace. Using E^ to approximate the slow subspace, n — n^ orthogonal- 
ity conditions are constructed and solved to compute points on the slow man ifold, 

shows 



1.2 



where n* is the dimension of i?" and the corresponding slow manifold. Fig. 
an example of a 2D slow manifold in M.^ and the relevant geometric objects. The 
spectrum shown is consistent with this geometry and would be constructed from the 
eigenvalues of D{{x) in the ILDM method. Kaper and Kaper ^T\ have analyzed the 
application of the ILDM method to a two-timescale system in standard singularly 



perturbed form ( 1.1 1 and shown that the error in determining the slow manifold is of 
order e^ and increases proportionally with the curvature of the slow manifold. The 
eigenvalues and eigenvectors of Di{x.) have also been used to express properties of 
flnite-time stable and unstable manifolds, under the assumption that the eigenvectors 
change sufficiently slowly with x along system trajectories [TB]. Eigenvectors are in- 
deed simpler to compute, than the Lyapunov vectors we will use, and should be used 
when they approximate the directions of interest to sufficient accuracy. Our focus 
however is on what to do when this is not the case (and for that matter, how to know 
when it is the case). 

Because the actual slow and fast subspaces E^{x.) and E^{x.) are invariant under 
the linear flow, it follows that in tangent vector coordinates adapted to these sub- 
spaces, the linear variational system must have a block diagonal structure such that 
the slow and fast dynamics are uncoupled. In general, the decoupling is not achieved 
using tangent vector coordinates corresponding to an x-dependent eigenvector basis. 
The computational singular perturbation (CSP) method [251 I26j includes an iterative 
procedure that adjusts the eigenvectors of D{{x.) to basis vectors that successively 
reduce the coupling between the slow and fast subsystems. In [3IT|, it was noted that 
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Fig. 1.2. Geometry of a two-timescale system with a 3D state space and a 2D attracting slow 
manifold. 



gap 



gap 




Fig. 1.3. Geometry of a two-timescale system with a 3D state space and a ID hyperbolic slow 
manifold. 



the basis vector refinement procedure in the CSP method is essentially a Lyapunov 
transformation used previously for subsystem decoupling in linear time- varying (LTV) 
systems [24] . Zagaris et al. ^44j have analyzed the accuracy of the CSP method ap- 
plied to a two-timescale system in standard singularly perturbed form (1.1) and found 



that the error in determining the slow manifold is of order e'^+^ after q applications of 
the CSP basis vector refinement algorithm |44j . The refinement requires accurately 
computing Lie derivatives of basis vectors |1D] . 

Both the ILDM and CSP methods rely on the eigenvalues of -Df (x) to determine 
the timescales and rely on eigenvectors to either determine (for ILDM) or initialize 
(for CSP) the tangent space structure. Rather than eigenvalues of -Df(x), existing 
theory |4l |22j for hyperbolic sets and normally hyperbolic invariant manifolds is based 
on Lyapunov exponents |28j - average rates of tangent vector length contraction and 
expansion along trajectories. Indeed the general invariant manifold theory developed 
by Fenichel [5] is based on Lyapunov type numbers, though eigen-analysis was appli- 
cable for his geometric singular perturbation theory because it is based on the local 
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structure for equilibrium points. In the case of a partially hyperbolic compact invari- 
ant set y C M" [TH], Lyapunov exponents are used to define an invariant splitting of 
the tangent space TxM" = E^''{x) ® E^{x) ® E^'^{x) into fast contracting, slow, and 
fast expanding subspaces at each x G y. If the slow distribution E'' is integrable, 
there is a corresponding fohation of the state-space. If a slow manifold S exists, it 
can be characterized as a leaf of the foliation that is invariant with respect to the 



flow, i.e., at each point on S, T^S — E^{ii.) and f(x) G £^'*(x). Fig. 1.3 illustrates 
a ID hyperbolic slow manifold in M.^ with the tangent space splitting shown at one 
point X S iS. The associated spectrum of Lyapunov exponents depicted in Fig. 1.3 is 



consistent with this splitting. Diagnosing and computing such geometric structure. 



encompassing both the attracting slow manifold (Fig. 1.2 ) and general hyperbolic slow 



manifold (Fig. 1.3 1 cases, is the goal of our work on two-timescale behavior. 



However, the geometric structure just described has been defined using asymptotic 
Lyapunov exponents. The Lyapunov exponents indicate average exponential rates 
over an infinite time interval; i.e., they are defined as limits (limit suprema in the 
most general case) when t — > ±oo. The state-space region is assumed to be compact 
and invariant under the flow. Under these assumptions, the infinite-time Lyapunov 
exponents converge, at least in the "lim sup" sense, and are metric independent. The 
tangent space splitting they define is invariant. But for many applications, the state- 
space region of interest is not invariant and the time span of interest is finite. The 
particular motivation for our work is to determine if disparate timescales are present 
in various fiight guidance problems, and if they are, to develop accurate reduced-order 
models to facilitate analysis and control design. The state-space region is the fiight 
envelope for the vehicle under study, and it is not invariant. The behavior outside 
this region may be different than that inside it, or the mathematical model may not 
even be valid outside the region, so we do not want the timescale information to 
be influenced by behavior outside the region. Hence the characteristic exponential 
rates of interest are averages over a finite-time interval, and we are led to the use of 
finite-time Lyapunov exponents (FTLEs) and the corresponding finite-time Lyapunov 
vectors (FTLVs) to characterize the tangent space structure of a two-timescale system. 
For timescale analysis applicable to the general case of a normally hyperbolic slow 
manifold, the FTLE and FTLV (FTLE/V for short) approach we develop seems to 
be most appropriate. If the slow manifold is either attracting or repelling, it can be 
discovered via simulation in forward or backward time respectively, although even in 
this case other approaches can be beneficial [5j [29l |39]. There are certainly other 
special cases, such as determining manifolds associated with equilibria or periodic 
orbits, where other approaches, for example ones based on eigen-analysis of Df, are 
applicable and may be more efficient. Chemical kinetics researchers have begun to use 
FTLE/ Vs [SS] , and eigen-analysis of the fundamental solution matrix $ for the linear 
variational equations over a finite time interval was considered in |39j . We note that 
the FTLE/Vs can be derived from eigen-analysis of <I>^<I>, reducing to eigen-analysis of 
the symmetric part of Df for infinitesimal propagation time, whereas eigen-analysis 
of <& reduces to eigen-analysis of Di for infinitesimal propagation time [7]. Eigen- 
analyses of Di, <&, and ^^^ for characterizing the fiow on an attractor were studied 
and compared in [[11 . Researchers in fiuid dynamics |42] have used FTLE/Vs to 
understand local fiow features; more recent work j7l 116} 138] has used the maximum 
FTLE as an indicator of finite-time stable and unstable manifolds. Another approach 
to stable and unstable manifolds in finite-time vector fields for fiuids has been taken 
in [37]. 
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We refer to the use of Lyapunov exponents and vectors, whether asymptotic or 
finite-time, to analyze the Hnear variational equations as Lyapunov Analysis. Pio- 
neering investigations of the properties of FTLE/ Vs can be found in [TTJ [T3] HU [57] . 
For the specific purpose of characterizing two-timescale behavior, a previous paper 
[31] focused on the properties of the FTLEs: their relationship to the infinite-time 
Lyapunov exponents and the kinematic eigenvalues, and their metric and coordinate 
dependence. Motivation from fiight mechanics for decomposing dynamics on the basis 
of fast and slow behavior, and the relationship of Lyapunov exponents and vectors to 
the geometry of singularly perturbed systems was described in [32] • In the present 
paper, we present, in terms of FTLE/Vs, a definition of, and a means of diagnosing, 
two-timescale behavior of a nonlinear, finite dimensional, time-invariant dynamical 
system on a non-invariant subset of M". This leads to a method of computing a slow 
manifold, when one exists. The efficacy of the method depends on the convergence 
rate of certain tangent subspaces defined by FTLVs as the averaging time increases. 
Previous convergence results [TT| [3T] are improved by characterizing the convergence 
in terms of the distance between the critical subspaces rather than in terms of the 
convergence of individual FTLVs. The scope of the present paper ends with the pro- 
cedure for identifying points on a slow manifold. A goal of future work is to obtain 
reduced-order models for the dynamics on the slow manifold and the boundary-layer 
dynamics. In the context of several application examples, presented at the end of 
the paper, a first step is taken towards determining effective numerical algorithms for 
implementing our methodology, directly comparing with other methods, and tackling 
progressively more complicated and higher-order problems. 

The authors recently learned of the work by Adrover et al. [H H] . Their approach 
to timescale analysis is based on Lyapunov exponents and corresponding tangent 
space directions and subspaces like our approach. There are several distinctions be- 
tween their work and ours. They view finite-time Lyapunov exponents and vectors 
as a means of approximating the asymptotic counterparts and draw from asymptotic 
theory to interpret their numerical results. They consider attracting slow manifolds 
only, provide numerical evidence for exponential convergence of the tangent space 
structure, and introduce two numerical methods for performing the computations. 
We characterize the usefulness of finite-time FTLE/Vs also when it is not appropri- 
ate to consider asymptotic limits and propose definitions of two-timescale behavior 
and a slow manifold suited to the finite-time setting. We derive an explicit expo- 
nential bound for the convergence rate of the tangent space structure that depends 
on the timescale gap. We consider, in addition to attracting slow manifolds, general 
hyperbolic slow manifolds, necessitating the intersection of forward and backward 
filtrations. We use existing numerical methods. We have included two application 
examples from [1 to illustrate our method and facilitate comparison with their work. 

The rest of the paper is organized as follows. In section 2, we define the dy- 
namical system to be considered and recall some definitions from geometry. Section 
3 covers Lyapunov analysis: first we define FTLE/Vs and describe their use for the 
identification of the tangent space structure; second the asymptotic theory of par- 
tially hyperbolic sets is described briefly; and finally, we address the convergence of 
the tangent space structure. In section 4 we define a finite two-timescale set and 
present the conditions satisfied by points on a finite-time slow manifold. The proce- 
dure for applying the approach is given in section 5. In section 6 several application 
examples are presented to demonstrate the use of our method and compare it with 
other methods. Conclusions are given in section 7. 
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2. Dynamical System Description and Relevant Geometry. The method 
we develop will be applied to a particular coordinate representation of a dynamical sys- 
tem. Denoting the vector of coordinates by x € M", 2 < n < oo, the x-representation 
of the dynamical system is 

x = f(x), (2.1) 

where the vector field f : M" -^ M" is a smooth function. The solution of ( |2.1[ ) for 
the initial condition x is denoted by x(t) = 0(t,x), where 0(i, •) : M" -^ M" is the 
i-dependent fiow associated with the vector field f and 0(0, x) = x. We assume that 
(p is complete on M" for simplicity in stating some of the results, but the methodology 
developed will only be applied on a subset of the state space and the properties of the 
flow outside this subset are irrelevant. 



The linearized dynamics associated with (2.1) are 



V = Df (x)v. (2.2) 

We will analyze the linearized dynamics to characterize the timescales in the nonlinear 
dynamics. An initial point (x,v) is mapped in time t to the point (x(/;),v(t)) = 
(0(i,x), <i>(i,x)v) where $ is the fundamental matrix for the linearized dynamics, 
defined such that <i>(0,x) = /, the nx n identity matrix. With this initial condition, 
we refer to the fundamental matrix as the transition matrix. Geometrically, for a pair 
(x,v), we view v as taking values in the tangent space at x. The coordinates for v 
correspond to a tangent space frame whose axes are parallel to those of the x-frame, 
but with origin at x. The tangent space at a point x e E" is denoted by TxM". 
The tangent bundle TM" is the union of the tangent spaces over M" and (x, v) is a 
point in the tangent bundle, with v the tangent vector and x the base point. We 
need the interpretation (x, v) G TM", because the analysis of the linearized dynamics 
will define a subspace decomposition in the tangent space and the orientation of the 
subspaces will vary with the base point x. 

We adopt the Euclidean metric for W\ Consistent with the Euclidean metric, we 
use the Euclidean norm to define the length of a tangent vector, i.e., for v S T^M", 
the length is ||v|| — (v, v)^/^ and (•, •) is the standard inner product. We will also use 
the distance between equidimensional subspaces 5*1 and 5*2 of M" given by [T5] 

dist(5i,52)-||Pl-P2||2, (2.3) 

where Pi and P2 are orthogonal projections onto Si and S2, respectively, and || ■ II2 is 
the induced 2-norm. The distance has a value in the interval [0,1]. The largest 
principal angle 9 G [0,7r/2] between the equidimensional subspaces is defined by 
sin6' = dist(S'i,S'2). 

Let Wi, . . . , Wfc, k < n, denote vector fields defined on M" that vary smoothly 
with x and have the property that at each x e M", the vectors Wi(x), . . . ,Wfe(x) are 
linearly independent in T^M.". Then at each x, A(x) = span{wi(x), . . . , Wfe(x)} is a 
fc-dimensional subspace. If fc = n, then A(x) — TxM" and for each x the set of vectors 
provides a basis for T-^M.". li k < n, then A(x) is a linear subspace of TxM" and A 
is called a distribution on M". A distribution is ^-invariant, if for any x G M" and 
V G A(x), the property <l>(i,x)v G A((/)(i,x)) holds. Distributions Ai,...,Am allow 
a splitting of the tangent space, if T^^:^^^ ~ ^i(^) ffi ■ • • ® Am(x), where ® denotes 
the direct sum of linear subspaces. If each distribution in the splitting is <l>-invariant, 
then the splitting is called a ^-invariant splitting. 
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A smooth submanifold A4 C M" of dimension m < n is 0-invariant, if for any 
X G A^, 0(i,x) G A4 for all t. An equivalent requirement for invariance is that 
f (x) G T-x^Ai for all X G A^. We mention three ways of representing such a manifold. 

1. Algebraic constraints: At least locally, A^ = {x G K" : /ii(x) = ••• = 
hn-rni'^) — 0} whcrc hi, i — I, . . . ,n — m are independent constraints and 
smooth functions of x. Given the invariance of A4, for all x G A4 the con- 
straint functions satisfy Lf/ii(x) = (^(x),f(x)) = 0, i = 1, ...,n — m where 
Lfhi denotes the Lie derivative of hi in the direction f . 

2. Graph of a function: At least locally, the coordinates of x can be sepa- 
rated into a vector yiindep of m independent variables and a vector x^ep of 
n — m dependent variables, and there exists a function 7 : M™ -^ M"^™ 
such that Al = {x G K" : x^ep — 7(xmdep)}- Given the invariance of M., 
the function 7 should satisfy fdep(xi„dep,7) = Q^^^^indepi^indep,!) where 

^dep — ^depxp^indepj^dep) f^HO ^indep — ^indepy^indepT^dep) c^re uennefl Consis- 
tently with X = f(x). 

3. Tangent space splitting and invariance-based orthogonality conditions: In this 
case, we assume that any other smooth invariant manifold containing M 
coincides with A4. For each x G M", suppose there is a splitting of the 
tangent space TxM" — r(x) (r(x))^ into m and n — m dimensional sub- 
spaces, where (r(x))-'- is the orthogonal complement to r(x), the distribu- 
tion r is ^-invariant, and for each x G A4, T^Ai — r(x). Given a basis 
{wi(x), . . . ,w„_m(x)} for (r(x))-'-, the manifold can be defined implicitly 

by 

7W = {xgM" : (wi(x),f(x)) =0, i = l,...,n-m}, (2.4) 

where (•, •) denotes the standard inner product. 



The orthogonality conditions for f in (2.4 1 can be viewed as partial-equilibrium 



conditions, partial in the sense that the vector field f need only be zero in certain 
directions. If one has constraint functions hi, i = I, . . . ,n — m for representation 
1, then (r(x))-'^ — span{ ^^ (x) , . . . , — g^-^(x)}. It may be easier however to find 
a basis for (r(x))^ directly without first finding constraint functions. Not every 
basis of (r(x))-'- can be related to a set of constraint functions. Determining the 
scalar constraint functions in representation 1 and the vector-valued function 7 in 
representation 2 requires the solution of partial differential equations and posing these 
equations requires a priori knowledge about the manifold, e.g., its dimension and in 
the case of representation 2, its orientation. 

The approach developed in this paper leads to a characterization of points on a 
slow invariant manifold as in representation 3. The ILDM approach does this also. 
The distinction is that the ILDM approximates the splitting r(x) (r(x))-'- based on 
the eigen-analysis of I?f(x), whereas our approach is based on finite-time Lyapunov 
analysis. 

3. Lyapunov Analysis. In this section we present the methodology for charac- 



terizing the linearized dynamics (2.2 1, along trajectories of the nonlinear system (2.1 1, 
that will enable the definition and diagnosis of two-timescale behavior. We refer to 
this methodology as Lyapunov analysis and think of it as serving the purposes for 
LTV dynamics that eigen-analysis serves for LTI dynamics. In the first subsection, 
we present a finite-time version of Lyapunov analysis, modeled after the asymptotic 
version described in Barreira and Pesin j4j and Katok and Hasselblatt |22j. We use 
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some of the notation and style of those books. In subsection 3.2, a brief account is 
given of how asymptotic Lyapunov exponents are used to define an (asymptotic) two- 
timescale set. In the final subsection, the convergence rate of a Lyapunov subspace is 
characterized, setting the stage for the finite-time approach presented in the remain- 
ing sections. Computational methods for Lyapunov analysis are considered briefiy in 
subsection 5.3; see [6] and references therein for the state of the art. 

3.1. Finite-Time Lyapunov Exponents/ Vectors and Tangent Space Struc- 
ture. A vector v € T:xM''' , propagated for T units of time along the trajectory 4>{t, x), 
evolves to the vector <i>(T, x)v in the tangent space T^i^t,x)^"- The ratio of the 
Euclidean lengths of an initial non-zero vector and its corresponding final vector, 
cr(T,x,v) = |j<I>(T,x)v||/||v||, is a multiplier that characterizes the net expansion 
(growth) if ct(T, x,v) > 1, or contraction if cr(T, x,v) < 1, of the vector over the 
time interval [0,r]. We distinguish variables associated with forward-time propaga- 
tion and backward-time propagation using the superscripts "+" and "^" respectively. 
The propagation time T, also referred to as the averaging time, is always taken to be 
positive whether forward or backward. The forward and backward FTLEs are given 
by 

M+(T,x,v) ^ ^lna+(r,x,v) = llnMJ>M, 

M-(T,x,v) = ^lna-(T,x,v) = ilnM^I|£M, ^^'^ 

for propagation time T. For v = 0, define /i+(T, x, 0) = /x~(r, x,0) = — oo. A 
Lyapunov exponent allows the corresponding multiplier to be interpreted as an aver- 
age exponential rate, i.e., (t(T, x,v) — exp[/i(r, x, v)r]; the average is over the time 
interval [0,r]. 

Discrete forward and backward Lyapunov spectra, for each (T, x), can be defined 
as follows. Define 1^{T, x), i — 1, . . . , n, to be the orthonormal basis of TxM" with the 
minimum sum of exponents, i.e., the minimum value of Y.f^iiJ,J{T,x,lf{T,x.)) over 
all orthonormal bases [SI. The forward Lyapunov spectrum is the set of exponents 
corresponding to the minimizing solution, namely, {fi^{T, x),i — 1, . . . ,n}. The Lya- 
punov spectrum is unique, though the minimizing basis is not in general. One way 
[^I31| to obtain a minimizing basis and the forward Lyapunov spectrum is to compute 
the singular value decomposition (SVD) of $(r,x) = N+{T,x)Y.+ {T,x)L+{T,x)'^ , 
where S"'"(r, x) = diag{ai{T,x), . . . ,a^{T,x)) contains the singular values, all pos- 
itive and ordered such that a^{T,x.) < CT^(r, x) < • • • < cr,^(T, x), and to compute 
the Lyapunov exponents as ^f{T,x.) = (l/T) ln(T^(T, x), i — l,...,n. The col- 
umn vectors of the matrix L+ (T, x) are the minimizing orthonormal basis vectors 
1^(T, x),j = l,...,n for TxM", and the column vectors of the orthogonal matrix 
Af+(r, x) are denoted n^(T, x),« = l,...,n. Rearranging the SVD of <i>(r, x), we 
can write $(r,x)l+(T,x) = exp[/i+(r,x)T]n+(T,x) which indicates that n+(r,x) G 
^0(T,x)II^"- Geometrically, the unit n- sphere centered at the origin in T-xM" propagates 
under the linearized dynamics to an n-dimensional ellipsoid in T^(T,xi)^""j the princi- 
pal semi-axes of the ellipsoid are exp[/i^(r, x)T]n^(r, x), i = 1,. . . ,n and the unit 
vectors in T^R"' that evolve to these vectors are respectively 1^(T, ii.),i — 1, . . . , n. 

Similarly, the backward Lyapunov spectrum consists of the exponents for the unit 
vectors in T-xM" that map to principal axes of an n-ellipsoid in r0(_7- x)!^"- The back- 
ward exponents can be obtained from the singular value decomposition $(— T, x) = 
7V-(T,x)S~(r,x)L^(r,x)'^ by M,^(r,x) = (l/r)ln(T^"(r,x), i = l,...,n. Assume 
the ordering on the diagonal of E^(r, x) is such that <Tj^(r, x) > ••• > (t,^(T, x). 
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The column vectors of the orthogonal matrix L (T, x) are denoted by \ (T, x), 
i = 1, . . . , n. For the column vectors of L^ (T,Ti.) and the orthogonal matrix N^ (T, x), 
we have l,"(T,x) e TxM" whereas n~(T,x) e r0(-T,x)IR"- 

In summary, a unit n-sphere in T^M." is propagated T units of time forward to 
an n-ellipsoid in T^[t,x)^^ and backward to another n-ellipsoid in Tij,(_t,x)^"' ■ In 
TxM", the 1^(T, x) vectors propagate to the principal axes of the forward ellipsoid, 
whereas the 17 (T, x) vectors propagate to the principal axes of the backward ellipsoid. 

The l+(T,x) and the 1^(7, x) vectors, for 



See Figure 3.1 
1 



for the case of n = 2. 
I = i, . . . , n, referred to as forward and backward FTLVs, respectively, will be used to 
define subspaces in TxM" associated with different exponential rates. Methods based 
on QR decomposition provide alternatives to computing FTLE/Vs [B]; see section 5.3. 




Fig. 3.1. Trajectory of nonlinear system and associated tangent spaces, illustrating the role of 
the Lyapunov exponents and vectors in the forward and backward propagation of a sphere of tangent 
vectors. Blue objects correspond to forward propagation, and green objects correspond to backward 
propagation. The arguments (T, x) of the FTLE/Vs have been suppressed. 



Definition 3.1. [N on- Degenerate Lyapunov Spectra] The forward and backward 
Lyapunov spectra are non-degenerate for particular arguments (T, x), if there are n 
distinct forward FTLEs and n distinct backward FTLEs, respectively. 

Assumption 3.1 For allT and^i. under consideration, the forward and backward 
FTLE spectra are each non-degenerate. 

This assumption simplifies the presentation. It will be modified for the subspace 
convergence proof presented later. We note however that distinctness is also related 
to integral separation and the stability of the Lyapunov exponents with respect to 
perturbations in the linearized system matrix, Df{x) [^. In some of the application 
examples presented in section 6, there are degeneracies; however, these degeneracies 
occur in an initial "transient" phase that is short relative to the time interval under 
consideration. Modifying the assumption to hold for T >To, for an appropriate value 
of To, is sufficient for applying the methodology we develop. 

The following subspaces, for i — l,...,n, can be defined by the orthonormal 
FTLVs 






— span{lf{T,ii.), 
= span{l7(T,x), 



,l+(T^,x)}, 



(3.2) 



and will be referred to as finite-time Lyapunov subspaces. For any i e {1,2, . . . ,n}. 
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//+(r, x,v) < ij'^{T,x) for any v G £^{T,x). However, for finite T, there also exist 
vectors v e T^M" \ C+{T,x) for which ^+(T,x,v) < nt{T,x). Although stated only 
for the forward-time case, analogous properties hold for the backward-time exponents 
and subspaces. 

If a collection of r < n linear subspaces of TxM" can be ordered such that Ai (x) C 
A2(x) C • • • C Ar(x) = TxM" with all inclusions strict, then this collection of nested 
subspaces defines a filtration of TxK". The nested sequences of subspaces 

{0} = £0 C C+{T, x) C /:+(T, x) C • • • C C+{T, x) = TxM", (3.3) 



TxM" - C^{T,x) D £^(T,x) D • • • D C-{T,x) D £"+1 - {0}, (3.4) 

are forward and backward filtrations [H [52] of T^xK^ . 

We need both forward and backward filtrations, because their intersections are of 
particular interest, as motivated by the following. Consider a two-dimensional non- 
linear system with an equilibrium point x^ . Assume the linearized dynamics at Xg are 
characterized by distinct eigenvalues A/ and A^, with A/ < A^ < 0, and correspond- 
ing eigenvectors e/ and eg. As T ^ 00, the FTLEs approach the eigenvalues, i.e., 
ix^ -^ Xf and /xj ^ As, and the first Lyapunov vector approaches the corresponding 
eigenvector Ij^ -^ ef. The second Lyapunov vector Ij approaches e^ the vector 
perpendicular to e/. The subspace £j'"(T, Xg) thus approaches E^{xe) = span{ef}, 
the eigenspace for A/ as T — + 00, whereas £j(r, Xe) — Tx^M^ for any T. It would be 
desirable instead to obtain the invariant splitting Tx^M^ = i?^(xe) © £'*(xe) where 
E"{xe) = span{es}. However, asymptotically all the vectors not in Ci have the Lya- 
punov exponent /i2 = ^s', thus the Lyapunov exponents for forward-time propagation 
do not distinguish £'*. The way to obtain E'' is by repeating the same analysis for 
backward-time propagation; in this case, the situation is reversed: asymptotically 
1^ -^ eg and E'^ can be distinguished, whereas E^ cannot J22ll43j. 

Because the Lyapunov exponent and vector information concerns how the lengths 
of vectors change, this information is in general dependent on how the length of 
a vector is measured. The FTLEs and corresponding tangent space geometry are 
invariant with respect to coordinate transformations provided the representation of 
the metric is transformed consistently. In the asymptotic case the Lyapunov exponents 
are also metric invariant, but this is not true for the finite-time case. This issue 
was addressed in |31j . In the present paper, we use the Euclidean metric exclusively, 
though any Riemmanian metric could be used. If two-timescale behavior is not present 
in the original metric under consideration, there may be another metric for which there 
is two-timescale behavior, as noted by Greene and Kim [14]. We are not addressing 
this opportunity directly, although one can apply our method with different metrics. 

3.2. Asymptotic Lyapunov Exponents and Two-Timescale Set. We draw 
from [18 to present a brief account of the asymptotic theory, covering only those 
definitions and results that serve to motivate and support our definitions and results 
for the finite-time case. 

A compact invariant set y C M" is a uniform two-timescale set (called uniform 
partially hyperbolic in '^E\), if there exists an invariant splitting at each x G y 

T^R" = Ef''{x)®E^x)®Ef''{x), (3.5) 
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and constants /i'', ^J, and C, with < /i" < /i-^ and C > 0, such that Vi > 

ve£;^=(x) => ||$(t,x)v|| < Ce-^^'*||v||, 

ve£;"(x) ^ C-ie-^°*||v|| < ||$(t,x)v|| < Ce^'*||v||, 

veEf^ix) ^ ||$(-t,x)v|| < Ce-^^'*|jv||. 

Consider a compact, invariant set y C M". When the infinite-time hmits (T -^ 



oo) of the exponents in (3.1) exist at x G 3^ for all v g T^R", they are denoted 
by /i~''(x,v) and /i^(x,v) and the system is said to be, respectively, forward regular 
and backward regular at x. There are at most n distinct exponents for the vectors 
in Ta;M"'\{0}. Consistent with our assumption for the finite-time case, we assume 
that there are n distinct exponents, denoted fi'[{x), i = 1, . . . ,n for forward-time and 
fi^ (x), i — 1, . . . , n for backward-time, with the forward exponents in ascending order 
and the backward exponents in descending order. Lyapunov subspaces are defined 
by C+{x) = {v e r,E" : /x+(x,v) < fi+{x)} and CHx) = {v e T,M" : ^l~Jx,^r) < 



/i^(x)}. Forward and backward filtrations are defined as in (3.3 1 and (3.4 1. The 
system is Lyapunov regular [3] at x if (i) it is forward and backward regular at x, 
(ii) /i^(x) = — /j,~(x), i = 1, . . . , n, (iii) the forward and backward filtrations have the 
same dimensions, (iv) there exists a splitting Txy = Ei{x)(B- ■ •0£'n(x) into invariant 
subbundles such that Cf{x) = Ei{x) (B ■ ■ ■ (B E^{x) and C~{x) = Ei{x) ® ■ ■ ■ & -E„(x), 
i = l,...,n, and (v) for any v G E,{x) \ {0}, limt^±oo(lA) In ||$(t,x)v|j = M+(x). 
The invariant splitting described in parts (iv) and (v) is referred to as Oseledet's 
decomposition . 

Next we describe how the asymptotic Lyapunov exponents can be used to char- 
acterize a two-timescale set. For the purpose of motivating the finite-time the- 



ory presented in the next section, we assume the system (2.11 is Lyapunov reg- 
ular at all the points of a compact, invariant set y, e.g., a periodic orbit. Fur- 
ther, assume that at each x G 3^, there are n-'"'^ large negative exponents, n* small 
in absolute value exponents, and n^'^ large positive exponents, with n-^'^ + n'^ + 
n-^^ — n. That is, there is a splitting of the forward Lyapunov spectrum of the 
form Sp+{x) ^ Spf%x) U Sp^{x) U Spf'^ix) where Spf'^ix) = {/i+(x), . . . , M+/c(x)}, 
Sp'ix) = {^+,^^(x),...,^+^,_^^,(x)}, and Spf'^ix) = {^+.,_^^,_^^(x), . . . ,/i+(x)} 
and constants < ii" < ^J defined by 

fif = min{-max/i+^,(x),min^+^,_^„,_^^(x)}, 
M^ = max{-min/i+^,^^(x),max/i+^,^,^,(x)}, 

where the maxima and minima are taken over the set y . If 3^ is a periodic orbit, 
the exponents do not depend on x and there is a zero exponent corresponding to the 
direction tangent to the orbit, but we will not restrict our discussion to the particular 
case of a periodic orbit. Using the forward and backward filtrations one can construct 



the invariant splitting (3.5 1 by 



E%^) - 4.+„.(x)n£;,,^^(x), (3.7) 



With fj.f and ^i" as defined in ^^, C = 1, and E^", E" , and E^" as defined in ( |3/7| , 
y is an asymptotic uniform two-timescale set. 

Although Lyapunov vectors were not used to define the subspaces in the asymp- 
totic case, they can be defined as follows and could be used to define the Lyapunov 
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subspaces. Let {lf{x.),i — 1, . . . ,n} denote an orthonormal basis for T^W^ such that 
{1^ (x) , j = 1 , . . . , i} is a basis for £+ (x) for every 1 < i < n. Let {1^ (x) , i = 1 , . . . , n} 
denote an orthonormal basis for T^M" such that {l~(x), j = i, . . . ,n} is a basis for 
£~(x) for every 1 < i < n. When there are n distinct Lyapunov exponents as we 
are assuming, it follows that these bases are unique up to multiplication of individual 
vectors by ±L These are clearly the orthonormal bases that minimize the sum of the 
asymptotic exponents over the set of orthonormal bases, and hence the basis vectors 
are the asymptotic counterparts of the FTLVs. 

Thus in the asymptotic setting either Lyapunov exponents or vectors (as just 
defined) can serve to define the Lyapunov subspaces and tangent space sphtting, and 
the results are equivalent. In contrast, the FTLEs and FTLVs define different tangent 
space structure. If one defines the «*'' forward finite-time Lyapunov subspace at x as 
C^{T,-x) = {v G Ta;K" : ^+(T, x, v) < ^f{T,ic)}, one gets not a subspace, but a cone 
centered on the FTLV-defined Lyapunov subspace £f{T,x). To see this, consider 
the tangent vector v = u + aw in T^M", with u € CJ{T,yi), w e (£+)-L(T,x), 
and a a scalar constant. For a given T, there exists a nonzero value of a small 
enough that v will belong to C^{T, x), whereas it does not belong to Cf{T, x). Under 
certain conditions, as T increases, £f{T,x) converges to its asymptotic value >C^(x) 
(as characterized in the next subsection) and Cf{T, x) converges to Cf{T, x) and thus 
to Cf{x) as well. 

Because the FTLV-defined Lyapunov subspace convergence is exponential in T 
(see next subsection), while the Lyapunov exponent convergence is much slower, per- 
haps proportional to 1/T [TT|, in the finite-time setting, we define the Lyapunov 
subspaces in terms of the FTLVs. 

3.3. Exponential Lyapunov Subspace Convergence. In this subsection, we 
relate the finite-time tangent space structure introduced in section 3.1 and the asymp- 
totic tangent space structure described in section 3.2. Theorem |3.5| below gives the 
exponential rate at which the finite-time Lyapunov subspaces introduced in section 
3.1 and expressed in terms of the FTLVs evolve with increasing T toward their asymp- 
totic limits. Most of the main ideas in Theorem |3.5| and its proof can be found in 
[llj . The new element here is that convergence of a particular Lyapunov subspace 
is addressed explicitly, rather than the convergence of individual Lyapunov vectors. 
It is this specific convergence rate property on which the methodology described in 
the following section rests. Before presenting the theorem, a couple definitions and a 
proposition are needed. 

The following proposition provides a formula for computing the distance between 
the subspaces £^(ri,x) and £J"(T2,x) in T^W^ for any value of j in the index set 
{1,2,. ..,n}. 

Proposition 3.2. Let L^{T,x.) denote the matrix whose columns are the Lya- 
punov vectors l^(r, x), i — 1,. . . ,j, and L'J,{T,x.) denote the matrix whose columns 
are the Lyapunov vectors lf{T,x), i = j + 1, . . . ,n. Then the distance between the 
subspaces £^{Ti,x) anrf £^(T'2,x) is 

dist{C+iT,,x),C+{T2,x)) =|li+(ri,x)^L+(r2,x)|l2 

= ||L+(r2,x)^L+(ri,x)|J2. ^^■^> 



Proof: Proposition 3.2 is a special case of Theorem 2.6.1 in [T^], page 76, and the 
facts that the columns of L'^{T,x) provide an orthogonal basis for £^(T, x) and the 
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columns of LJ', (T, x) are mutually orthogonal to the columns of L'J (T, x) . D 

Definition 3.3. Ul^ The Lyapunov spectrum is strongly non-degenerate at a 
point X, if there exists positive constants To and S such that the spectral gap between 
each neighboring pair of forward FTLEs, fj,f_^j^{T, x) — nj {T, x), j = 1, . . . , n — 1, is 
greater than 5 for all T > To and likewise for the backward exponents. 

To consider the convergence of a Lyapunov subspace C~^ (T, x) with T, we focus 
on a particular spectral gap and define it more precisely. 

Definition 3.4. [Relative Spectral Gap] For a specified To > 0, the relative 
spectral gap A/i (x) between neighboring forward FTLEs iJ,^{T,x.) and //^, j(r, x), 

for a particular j £ {1,2, ... ,n — 1}, is A/i~''(x) = miT>To{lJ-t+i{T,x) — fi'l{T,x)). 
The relative spectral gap A/j,~(x) between neighboring backward FTLEs ^[^^^{T^yi) 
and /ij7 {T, x) is similarly defined. 

Theorem 3.5. Consider the dynamical system (2.1) on a compact invariant 
subset y of the state-space M". At a Lyapunov regular point x G 3^ for which there 
exists To > and 6 > such that the Lyapunov spectrum is strongly non-degenerate 
for T > To and for which there is a nonzero relative spectral gap A/i (x) for a specific 

value of j , the subspace £~^{T,'k) approaches the fixed subspace Cf{x), defined in 
section 3.2 in terms of the asymptotic Lyapunov exponent /i^(x), at an exponential 



rate characterized, for every sufficiently small AT > 0, by 

dist{C+{T,x),C+iT + AT,x)) < i^e^M'")-^, 



(3.9) 



for all T > To, where K > is AT dependent but T independent. Similarly, as 
T increases, the subspace C^{T,yi) approaches the fixed subspace £j^(x) at a rate 
proportional to ea:p(— A/ijr(x) • T) where A/i^(x) = infT>T„(/ifc_i(r, x) — fi'f^{T,x)). 
Proof: Using (3.2) we have 



dist (/:+(r,x),/:+(T + AT,x)) = |iL+(T,x)^i+(r + at,x)|12 

[i++i(r + Ar,x) ... i+(r + AT,x) ] 






irm-)' 



(i+(r,x),i++i(r + AT,x)) 



(l+(T,x),l+(T + AT,x)) 



(l+(T,x),l++i(T+AT,x)) ... (l+(T,x),l+(T + AT,x)) , ^ 

(3.10) 
Borrowing a result from [llj . we have for T > to I'^'-order in the time increment 
AT 

l+(T + AT) = (l + cAT)l+(T) + AT ^ ^:..r ..L ..+ ".1.1 . (S-H) 



--lii^r. 



eit^t-t^t)T 



e(M+-M;t,)T' 



where A ~ T)f(x) is the system matrix of the linearized dynamics (2.2), n^ is a vector 



from the SVD of the transition matrix <i>(T, x) as defined in section 3.1, c is a constant 
that is inconsequential in the following developments and is thus left unspecified, the 
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X dependence has been suppressed, and all exponents and vectors in the summation 
on the right-hand-side are evaluated at (T, x). It follows that the inner products in 



(3.101 are 



HUT, X) , 4(T + AT, X)) ^ AT JI;:lllf_ +;,)<] ^ . (3.12) 

Because fee {1, . . . , j} and m e {j + 1, . . . ,n}, we have exp[(/i^(r, x) — ^+ (T, x))T] < 
exp[— A/Lt (x)r]. Let a — max2.£j;maxig{i 2,...,n} |'^i(^"'" + ^)|, the maximum eigen- 
value magnitude of A^ + A over the set 3^. And let a = exp(— 2A^"*"(x)ri) for some 
Ti>To. Then for T > Ti > we have 

|(l+(r,x),l+(T+ AT,x))| < ^e-^;W^. (3.13) 

Upper-bounding the 2-norm by the Frobenius norm and taking K = \J j(n ~ j)f3^, 
the bound in the theorem follows. This bound is conservative, due to the use of the 
Frobenius norm, but it shows the exponential rate of convergence. Using the bound 



(3.9), one can show that the sequence of iterates is Cauchy. Moreover this is true for 



every sufficiently small AT. Because the space of j-dimensional subspaces in T^W^, a 



Grassmannian, with the distance given in (3.2 1 as the metric, is complete, we conclude 



that £'^(T, x) approaches a fixed subspace. This subspace is Cf{'^) defined in section 



'J \ ^ / -ff ™ " - -i- -- -t-- - ~j 



3.2, because all vectors in it have exponents less than or equal to A't(x) and one can 
show that any vector not in the subspace must have a larger exponent. The proof for 
backward-time is similar. D 



The Theorem |3.5| hypothesis that the Lyapunov spectrum is strongly non-degenerate 
is necessary because the proof is based on the evolution of the individual Lyapunov 



vectors according to (3.111. We conjecture that the existence of the relative spec- 
tral gap is sufficient for the exponential subspace convergence, even if the rest of the 
spectrum has degeneracies. 

4. Finite-Time Two-Timescale Set and Slow Manifold - Theory. We 

define finite-time two-timescale behavior by first defining a finite-time uniform two- 
timescale set. A two-timescale set has a special tangent space structure, and allows 
us to formulate conditions that would be satisfied at points of a slow manifold. If 
a slow manifold exists, then the nonlinear system has two-timescale behavior of the 
boundary-layer type and there is an opportunity, though not pursued in this paper, 
for system decomposition and reduced-order analysis. 

We consider the timescale behavior of a system on a set X C M" which is in general 
not (/)-invariant. For the purpose of defining and diagnosing two-timescale behavior, X 
could be a point, a collection of isolated points, a segment of a trajectory, as examples, 
but in the search for a slow manifold, we assume no a priori knowledge and typically 
need to consider a bounded, connected open set of the state space. 

In the finite-time setting, the terms "fast" and "slow" are defined by qualitative 
properties of the dynamics, relative to a particular time duration T^, namely, "fast" 
refers to behavior that decays, either in forward or backward time, to a "negligible 
level" over Tc, whereas "slow" refers to behavior that does not. Although in each 
particular application, one needs to define fast and slow quantitatively, there is no 
generally appropriate definition; so we do not offer one. The bound [3 in the following 
definition is the means of quantitatively distinguishing fast from slow. A numerical 
value of j3 needs to be specified as appropriate for each application. We give some 
guidelines in section 5.1. 
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4.1. Finite-Time Two-Timescale Set. 



Definition 4.1. A set X C M", n > 2, is a uniform two-timescale set for (2.1), 
if there exist real numbers fi^ , jjl^ , To and Tc, with < /i* < /!•' and < To < T^, 
such that, on [To,Tc\ x X, there is a uniform splitting of the forward and backward 
Lyapunov spectra into fast contracting, slow, and fast expanding subsets, separated by 
gaps of size Afj, = fj,^ — fi"^ , where Afj.{Tc — To) > [3 and the positive constant (3 is a 
specified lower bound for two-timescale behavior. Specifically, the FTLEs satisfy the 
following properties for all (T, x) G [To,Tc] x X : 

1. M+/e(T,x) < -^lf, -Ai^ < /i+,^^(r,x), /i+,^„,(T,x) < ^^ and ^^ < 

2. -M;/c(r,x) < -iJif, -iJi^ < -m;,.+i(T,x), -m;,.^.„.(T,x) < /i^ and ^if < 

where n^'^, n^ and n-^"^ are constant positive integers, with n''^ + n^ + n'^ = n, that 
specify the number of exponents associated with fast contracting, slow, and fast ex- 
panding behaviors, respectively. Either n^'^ or n''^ is allowed to be zero, but not both. 
For n^'^ — 0, the conditions on iJ,'^f^,^^,^{T,x) and ii^f^ -^^{T,x) do not apply; 
similarly, for n^'^ = 0, the conditions on /x^j.^(T, x) and fj,~f^{T,x) do not apply. 

Properties 1 and 2 are illustrated in Figure [4T] where the bounds and forward and 
backward exponents are plotted on aligned different copies of the real line for clarity. 
The exponents for particular values of T and x are pictured, but Properties 1 and 2 
require this structure for all (T, x) E [To, T^y.X. The use of times up to T^ means that 
the averaging in computing the Lyapunov exponents and vectors involves trajectories 
which, though they begin in X , extend into the larger (unless X is (/)— invariant) set 

Xext = {y e M" : y = 0(i, x) for some (t, x) G [-T,, T,] x X}. (4.1) 
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Fig. 4.1. Spectra of forward and backward FTLEs illustrating the gaps (at the point x £ A" for 
the averaging time T). 



In Definition |4.1 1 the fast and slow behaviors are characterized by the exponent 
bounds /i-^ and /i'' respectively. To simplify the definition we have used the same mag- 
nitudes /i* and ^f to define both gaps; however the symmetry of the two gaps with 
respect to zero is not necessary. Properties 1 and 2 ensure that, uniformly in T and x, 
the forward Lyapunov spectrum can be divided into fast contracting, slow, and fast 
expanding subsets (Property 1) and the backward Lyapunov spectrum can be divided 
into fast contracting, slow, and fast expanding subsets (Property 2). Properties 1 
and 2 also ensure that common gaps in the forward and backward Lyapunov spectra 
not only exist, but also separate the spectrum in a dimensionally consistent manner. 
Although Properties 1 and 2 only apply to the exponents for points in X , they imply 
uniform timescale structure over X^xt, because the exponents represent averages over 
this larger set; in particular, the "kinematic eigenvalues" [HIST], whose averages pro- 
duce the exponents, must be similar over X^xt. To provides a grace period over which 
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the bounds on the exponents do not have to be satisfied, in order to accommodate 
initial transients which could otherwise violate the bounds. T^ is the characteristic 
maximum time over which the uniformity in the exponents holds and also is the time 
duration relative to which the terms fast and slow are relevant. One could use two 
difi^erent times; we have used Tc for the dual role to simplify the definition. The bound 
(3 also has a dual role: (i) it ensures that the fast motion decays to a negligible level, 
in either forward or backward time, over the duration T^, because the decay factor is 
exp(-^/(rc - To)) and is small if A^(rc - To) = {^If ~ Ai')(T; - To) > (i, and (ii) it 
ensures that the gap is large enough (once its value is specified) relative to T^ — To 
that the critical subspaces can be resolved, as clarified next. We note that if there are 
additional gaps in the spectrum, a multiple-timescale set could be defined similarly. 

Proposition 4.2. When X is a uniform two-timescale set, at each y: G X, the 
subspaces 

Ef%T,^) = /:+,.(T,x), 

E^{T,^) =. £+,.^„.(T,x)n£;,,^,(T,x), (4.2) 

i?/'=(r,x) ^ £;,.^„.^^(T,x), 

have, for all T £ [To,Tc\, the properties: 

V e Ef''{T,Ti)\{0} ^ ||$(r,x)v|| <e-^^^||v||, 

ve-E;'*(r,x)\{o} ^ ||$(r,x)v|| <e^'^||v| 



(4.3) 
and in addition, the subspaces E-^'^{T,x), E^{T,x.) and Ef'^{T,x.) approach, in the 



V e E'{T, x)\{0} ^ \M-T, x)v|| < e'^"^||v| 

V G Ef^{T, x)\{0} ^ ||$(-T, x)v|| < e-'^'^lh 



sense of (3.9), fixed subspaces, with increasing T, at least at a rate proportional to 
exp{—Aii ■ T), where A/i ~ ji^ — /i*. 

Proof: The four properties in (14. 3| follow from the definitions of the subspaces 



£+,.(r,x), £-,,^^(T,x), /:+,^„,(r,x), and £;,,^„,^^(T,x) (see (3.2)) and of a 



uniform two-timescale set (see Definition 4.1). Given the properties oi a uniform 



of Theorem 



3.5 



two-timescal e se t in Definition |4.1| an argument similar to that used in the proof 
can be used to show that /:+j,,(r,x), £;;j,^^(r,x), £+j,_^^,(T,x), 



and C f^, 3 .(T, x) approach with increasing T fixed subspaces at least at a rate 



proportional to e^^^^ . Because of the relationships in (4.2), it follows that E^'^{T,x), 
E'^{T,x.) and E^'^{T,x.) approach fixed subspaces at the same rate. Due to the finite- 
time constraint, this is not to say that convergence is achieved. D 

Proposition |4.2| says that there is tangent space structure associated with the 



slow and fast exponential rates. However the decay/growth bounds (4.3) hold for the 
subspaces £^^^(r,x), £;"(T,x) and Ef^T^yi) for all T G [To,rj. Which value of T 
should we use, i.e., which subspace structure is most appropriate? Proposition |4.2| 



states also that the subspaces are approaching fixed subspaces as T increases. If the 



hypotheses of Theorem 3.5 were applicable (if A" is a subset of a compact invariant set, 
etc.) and the T ^ oo limits could be computed, then we could compute E^'^{T,x), 
E''{T, x) and E^'^{T, x) at each point of X for arbitrarily large averaging times T and 
these subspaces would converge to <l>-invariant subspaces that depend only on x |4]. 
Limited to T G [To, Tc] we should use T = Tc to obtain subspaces that not only have 



the decay/growth bounds (4.3 1 but also approximate invariant subspaces as closely as 
possible within the available averaging times. If, however, An{Tc — To) is larger than 
the prescribed /?, it is sufficient to use the value of T satisfying A/i(T — To) — (3. 
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4.2. Finite-Time Slow Manifold. X being a finite-time uniform two-timescale 
set establishes the potential for the existence of a slow manifold. To define a finite-time 
slow manifold, we now assume X is an open set of M". 

Definition 4.3. A finite-time slow manifold is a submanifold of X denoted S{T) 
such that f(x) G E''{T,yi) for all x G S{T). 

The set 

{x e A- : (f(x),w) = 0,Vw e [E'{T,x)]^}, (4.4) 

thus satisfies a necessary condition and constitutes a candidate finite-time slow man- 
ifold. The motion at each point in the set is slow, because f(x) S i?^(x), and thus 
there is no fast motion, the components of f(x) in i?^'^(x) and E^'^{x) being zero. 
That this set is a manifold has to be verified. If a finite-time slow manifold ex- 
ists, it will in general not be relatively invariant with respect to X , because T,j;S{T) 
and E''(T,x) at points of the set do not coincide, and thus the invariance condition 
f(x) e TxS{T) is not satisfied. However, if _E*(r, x) is close to the (hypothetical) 
^-invariant asymptotic limit £'*(x), we conjecture that S{T) will be close to the cor- 
responding (hypothetical) i/i- invariant slow manifold. The examples in section 6, in 
which the asymptotic limits are relevant and can be determined, support this conjec- 
ture. 

5. Finite-Time Two-Timescale Set and Slow Manifold - Procedure. If 

the goal is only to diagnose two-timescale behavior and determine the tangent space 
structure, then X can be any subset of M". For example one could take X to consist 
of a single equilibrium point (a fixed point of the vector field) , although eigen-analysis 
would be applicable and more efficient for this particular case. If one also wants to 
search for a slow manifold, then usually X is an open set, because it will be necessary 
to iteratively search for points that satisfy slow manifold conditions in a state space 
region of full dimension. As an example, consider a normally hyperbolic periodic 
orbit in M.^, for which the transverse motion is faster than the motion along the 
periodic orbit. Let X be the neighborhood of a segment of the periodic orbit. The 
invariant slow manifold in X is the segment of the periodic orbit; the finite-time 
slow manifold, for a given T^, would approximate the invariant slow manifold. In the 
application of the methodology we do not require that X is embedded in an invariant 
set; this hypothesis was only used to relate our approach to the asymptotic theory 
and characterize the Lyapunov subspace convergence. 

5.1. Diagnosing a Finite-Time Two-Timescale Set. FTLEs are computed 
for a grid of points on X to determine if A" is a two-timescale set according to Definition 
|4.1[ One needs to see a pattern as illustrated in Figure |4.1| uniformljQ in T and x 
and to verify that the spectral gap is sufficiently large relative to Tc- Regarding 
uniformity, the individual exponents can vary with T and x as long as there is a 
sufficiently large uniform gap. Thinking of the reciprocals of ^-^ and A/Lt as time 
constantaM a guideline is that the fast behavior should decay/grow over several time 
constants and the finite-time subspaces E^'^{Tc,x), E^{Tc,ii.), and E^'^{Tc,k) should 
converge over several time constants toward their hypothetical infinite-time limits. 
If decay/growth and convergence over m time constants is desired, then we require 



* Because the FTLEs are only examined on a grid and at a finite set of values of T, some 
experimentation witfi the grid (in T and x) is required to ensure that it is sufficiently fine. 

TFor an exponential function of time, e~*' '^, the time constant t > is the time t at which the 
function equals e~^. 
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fi* (Tc — To) > m and A/i(rc — To) > "t-; the latter requirement is the most demanding, 
and in applying Def. |4.1| we should set (3 = m. For a given value of f3, the smaller 
the gap is, the larger Tc must be. However, unless X is 0— invariant, the set X^xt (see 



(4.1)) grows with Tc and at some point the timescale behavior may not be uniform 
on this extended set. 

We note that the convergence of the subspaces can be checked directly by moni- 
toring the distance between the subspaces with increasing averaging time (illustrated 
in section 6). To fill in the timescale structure around the grid points, one can either 
compute the FTLE/Vs at additional points or use the equations for propagating the 
Lyapunov subspaces (actually the orthonormal basis given by the Lyapunov vectors) 
along a trajectory given in |13j . 

5.2. Computing the Finite-Time Slow Manifold. Provided that X satisfies 
Definition |4.1| we can take the next step which is to look for a slow manifold in 
X, where X is now assumed to be an open set in M". Within X, the points on a 
candidate finite-time slow manifold S(T) are defined implicitly by the orthogonality 



conditions in (4.4 1. Rather than use local eigenvectors to form an approximate basis 
for the orthogonal complement to i?" as in the ILDM method, we propose using the 
appropriate Lyapunov vectors to form the basis for {E'^)-^ as given by the following 
proposition. 

Proposition 5.1. The orthogonal complement of E''{T,x) can be represented as 

(i?^(T,x))^ - 5pan{lr(T,x), . . . , r,,(T, x), l+,.+„.+i(r,x), . . . ,l+(T,x)}. (5.1) 



Proof: From Proposition 4.2 we have defined the slow subspace E^{T, x) = ^C^j^ , ^ (T, x)n 



£^;^,_^^(r,x). Using an identity ^, we have {E%T,:>i))-^ = (/:,7^,^^(T,x))-L ® 
(£^j,^ (T, x))-'-. The proposition then follows from the facts: (£^j,^ (T, x))-'- = 
spL {lr(T,x),...,r,,(T,x)} and (/:+,.^„.(r,x))^ = span {l+,"„,^^(r,x), . . . , 
l+(r,x)}.D 

In order to obtain solutions of the algebraic equations, we designate n* compo- 
nents of X as independent variables, fix their values, and determine the values of the 
remaining n — n'* components, the dependent variables, that minimize 

fc 

n-' n 

J = ^(lr(T,x),/(x))2+ Y. (l,V+„=+i(r.x),/(x))2, (5.2) 

j=l nf^+n^ + l 

for the particular value of T that has been chosen. Ideally at points on a finite- 
time slow manifold, the minimum value of J would be zero, but we use a numerical 
iterative solution procedure that is stopped once J is below a specified tolerance. 
This is repeated for a grid on the independent variable space. The directions of the 
Lyapunov vectors indicate how to separate the coordinates of x into independent 
and dependent variables. The independent variables must be chosen such that their 
coordinate axes are not parallel to any directions in (E^) . Different independent 
variables might be required in different regions of X. Within X, there could be zero, 
one, or more than one slow manifold. Hence, for fixed values of the independent 
variables, there could be zero, one or several minima. The set of points satisfying 
the orthogonality conditions must be further examined to identify whether there is a 
finite-time slow manifold. 

This procedure for determining a slow manifold is much the same as the ILDM 
method, except for the important difference that FTLE/Vs are used instead of eigen- 
values and eigenvectors (in practice ILDM is often implemented using Schur vectors. 
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an orthogonal basis generated from the eigenvectors) . Determining the Lyapunov vec- 
tor basis for (-E^)^ requires more computation than determining the eigenvectors, but 
the potentially greater accuracy may be needed, if the timescale separation is mod- 
est or if the slow manifold has significant curvature, since these are the conditions 
[2T| for which the accuracy of the eigenvector-based ILDM is reduced. The improved 
accuracy is demonstrated in the following section. Whether or not improved accu- 
racy can indeed be achieved in a numerical implementation of our method on more 
complicated and/or higher-dimensional applications than those considered in the next 
section remains to be determined. 

5.3. Numerical Methods. The numerical methods used for the computations 
presented in the next section are described in this subsection. All the computations 
are done in the Matlab® environment. The numerical integration of the nonlinear 
state equations and the corresponding linear variational equations is performed with 
the 'ode45' integrator. 

The FTLEs and FTLVs associated with an initial state x are computed for an 
averaging time T either by SVD or QR factorization. Only the computation of the 
forward-time FTLE/Vs is described, since the computation of the backward-time 
FTLE/Vs is analogous. The first step of both methods is to integrate the nonlinear 
state equations from t = to t = T and save the values of 0(i, x) at the N equally 
spaced times At, 2At, ..., NAt, where NAt = T. 

In the SVD method, the transition matrix is computed and then the SVD is 
applied. The transition matrix is computed by integrating, simultaneously, the non- 
linear equations and the linear variational equations over each segment of the base 
space trajectory, with the state initialized with the saved value at the beginning of 
the segment and the transition matrix initialized with the identity matrix. Using the 
notation <i>^* = <i>(A<, (^((/c — 1)-At,x)) for k = l,2,...,Af, the transition matrix is con- 
structed from the transition matrices for the segments as <i>(T, x) = <i>^* • • • <i>^*(f>^*. 
The resulting transition matrix is then factored as <i>(r, x) = A^+I]+(L+)-^ using the 
'svd' command in Matlab®. Each FTLE, /ii(T, x), is obtained by ^^^(T, x) = ^ In af , 
where ai is the i^^ singular value of $, the positive square root of the i*'* diagonal 
element of S"'". If this procedure does not produce FTLEs in the ascending order we 
have assumed in our notation, the procedure is modified to achieve ascending order. 
The FTLVs 1^(T, x), i = 1, . . . , n are the column vectors of L+. 

For a given trajectory from x to 0(r, x), for a particular T, we have the op- 
tion of computing certain Lyapunov vectors at x and at (j){T, x) by forward or back- 
ward integration. Because <I>(-T, 0(T, x)) = <i>(r, x)^\ it follows that i+(r,x) = 
7V-(T,0(r,x)) and7V+(T,x) = L-(r,(^(r,x)). 

In the QR method, a segment by segment approach is also used [S]. For the 
fc*''-segment, after the transition matrix is computed as described in the previous 
paragraph, the Qk~i matrix associated with the state at the end of the previous 
segment is propagated by the transition matrix to the end of the fc*''-segment and the 
QkRk factorization of the resulting matrix is obtained, as summarized by 

<^t'Qk~i = QkRk (5.3) 

This sequence of operations for k = 1,. . . ,N must be initialized by prescribing Qo] 
typically the identity matrix is used [BJ [TD] . It then follows that 

$(T,x)Qo = Q(T,x)i? (5.4) 
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where Q{T,x.) = Qm and R = i?Ari?7v-i ■ • • ^2^1- Note that if we choose Qo — 
L+(r,x)), then Q(T,x) = iV+(T,x), or equivalently (5(T,x) ^ L-{T,(j){T,yi)), and 
i? = S+. For almost every Qo, Q{T, (t>{T, x)) will approach N+{T, (j){T, x)) and R will 
approach E"*" in the absence of numerical errors. We have found the QR method to be 
more numerically reliable. In the following section, both methods produced essentially 
identical results for the 2D and 3D examples. For the 5D example, we found that the 
QR method allowed longer averaging times without exhibiting numerical problems. 

Adrover et al. [TJ |5j have considered two alternative computational approaches. 
One is based on exterior algebra, that calculates, in addition to the FTLEs, the Lya- 
punov subspaces directly, rather than calculating vectors that span these subspaces 
as in the SVD and QR methods. The other propagates bases for the tangent and co- 
tangent spaces along trajectories, with periodic adjustments to keep the basis vectors 
from all aligning with the direction associated with the largest Lyapunov exponent. 
We do not have any computational experience with either method. 

The approximation to the slow invariant manifold is computed from the orthog- 
onality conditions. We fix the values of n^ components of the state vector, x, as 



the independent variables, and solve the unconstrained minimization problem (5.2 1 
to find the values of the remaining n — nf components using the 'fminunc' function 
in the Matlab® Optimization Toolbox. 

6. Application Examples. Several application examples are presented to demon- 
strate the use of the finite-time Lyapunov analysis (FTLA) method and compare it 
to the ILDM method. 

6.1. Two-Dimensional Linear Time- Varying Example. The purpose of 
this example is to distinguish clearly between the timescale information provided 
by Lyapunov analysis and the timescale information provided by eigen-analysis. Con- 
sider the LTI system 



w = Aw 



\f 
A., 



(6.1) 



where w = (wi,if2)'^ and the eigenvalues of A are real with A/ < Ag < 0. We 
introduce the coordinate transformation 



i?(f)w 



cose{t) - sin 6l(i) 
sin e[t) cos e{t) 



(6.2) 



where 9 = cut and lu is a, constant. In terms of x, the system is 

X = A{t)x = {RR^ + RAR^)x = R{R^ R + A)i?'^x. (6.3) 

The solution for the w-system for the initial condition w(ii) = Wi is w(i) = exp(A(i— 
ii))^! where exp(A(i — ti)) is the transition matrix. The corresponding transition 
matrix for the LTV x-system is 

$(i, h) = R{t)e^^'-'''>R^{h). (6.4) 

In the w-representation the behavior is composed of fast and slow exponentially 
contracting modes with fixed directions given by the eigenvectors of A. In the x- 
representation there are also fast and slow exponentially contracting modes, but the 
fast and slow directions are rotating. We will show that the Lyapunov exponents and 
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vectors identify the slow and fast exponential modes, whereas in general the eigenval- 
ues and eigenvectors for A(t) and $ do not. This example serves as an idealization of 
the linearized dynamics of a nonlinear system whose slow and fast directions rotate 
along a trajectory, as would be the case along a trajectory on a slow manifold with 
curvature. To be more consistent with the notation of the previous section, we use 
T to denote the propagation time and let t denote the time at which the Lyapunov 
vectors are computed. 

The transition matrix <f>(i + T, t) for the x- system has the SVD 

$(t + T,t)^ N+{t + T)T.+{t + T, t){L+(t)f ^ R{t + T) exp(AT)i?^(i). (6.5) 

The Lyapunov exponents are /i^ = A/ and ^2 = As and are independent of t and 
T for this example. The Lyapunov vectors, independent of T, are the columns of 
L+(i) = L~{t) = R{t). The SVD of $ thus identifies the exponential rates of the two 
modes in S"*" and the rotating directions of these modes in L+ and L^ . Specifically we 
have the fast subspace Ef{t) = £i{t) = span{l^{t)} and the slow subspace E''{t) = 
£^(i) — spari{l^(i)}, where the fast direction is lf{t) = [cos9{t) sm9{t)]'^ and the 
slow direction is 1^ (i) — [— sin 9{t) cos 9{t)]'^ . In this example the rotational motion 
is periodic and Floquet theory is applicable; but, the SVD would also characterize 
different and irregular rotations of the fast and slow directions. 

The eigenvalues of A{t) — {RR^ + RAR^) are |3j the same as the eigenvalues of 
the matrix 



(R^R + A) 



A 



/ 



A. 



(6.6) 



because the two matrices are related by a similarity transformation. The two eigen- 
values of A{t) are 



i(A, + A,.±((A.-A,)2-4c.^)^/^), 



(6.7) 



and are denoted by A+ and A_ based on which sign is taken. The corresponding 
eigenvectors of A{t) are 



v+ = Rit) 
V = R{t) 



—w 
X+-Xf 



A- - Xs 



(6.8) 



The eigenvalues oi ^{t + T, t) are 

^ (cos{LoT){e^i^ + e^=^) ± (cos\ujT){e^'^ + e^^^f - 4e(^/+^=)^ ) ' ) . (6.9) 



1/2 



The eigenvalues of A and $ depend on the rotation rate oj. For A, the eigenvalues 
are independent of t and T for this example. For small lo, the eigenvalues of A 
are good approximations of the Lyapunov exponents and indicate the slow and fast 
behavior. However, as oj increases from zero the eigenvalues approach each other, 
becoming equal when uj — (As — A/)/2, and then split as a complex conjugate pair. 
The eigenvalues of $ depend on coswT. For a given value of lj, at times for which 
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COS LoT — 0, the eigenvalues are pure imaginary; only for T such that cos loT — 1 
are the eigenvalues e'^^^ and e^"'^ . Hence the eigenvalues of A and $ are not always 
reliable indicators of the fast and slow exponential modes. The eigenvectors of A 
provide the directions associated with distinct exponential rates only if these directions 
change slowly enough. 

6.2. Davis-Skodje 2D System: Attracting Slow Manifold. Davis and 
Skodje (DS) |[5j introduced a two-dimensional nonlinear system for which the slow 
manifold can be determined analytically. The DS system is 



Xi 



X2 



-Xi 



I (7-l)a:i+7^i 
72^2 + (l+a;i)2 



(6.10) 



defined on the state space {(a;i,a;2) e 



ci > and X2 > 0} with constant 



7 > 1. The origin is a globally attracting equilibrium point, but more importantly in 
the present context, for sufficiently large 7, trajectories are first attracted on a faster 
timescale to the ID slow manifold 

S = {(xi, xa) e M^ : X2 - a;i/(l + xi)}, (6.11) 

and then follow S to the origin on a slower timescale. The two timescales are evident 



in the analytic solution for the flow associated with the vector field in (6.10) 



{t;xi,X2) 



xie 



{-^ - ifk) 



-It 



(6.12) 



Note that if the initial state is on the slow manifold, there is no fast timescale 



behavior because the coefficient of e '>'* in (6.121 is zero. The slow manifold S and 



several other trajectories are shown in Fig. |6.1| for 7 = 10. The time interval between 
asterisks on the trajectories is 0.1, illustrating faster motion off S than on S. From the 



analytical representation (6.11 1 for the slow manifold, we know that for any x e 5, 

T^S ^ span{[{l + xif if}. (6.13) 

The linearized dynamics of the DS system are given by v == Z?f (x)v, where 



Dt 



-1 

(7-l) + (7+l)a;i 
(1+2:1)=' 





-7 



(6.14) 



Given the presence of the equilibrium point, other approaches based on eigen- 



analysis at the equilibrium point are applicable: for example, integrating (6.10 1 back- 



ward from an initial state perturbed slightly from the origin in the direction of the 
eigenvector associated with the largest eigenvalue to compute S. However our purpose 
here is to demonstrate the methodology developed in this paper, methodology that 
does not require the presence of an equilibrium point. 

6.2.1. Finite-Time Lyapunov Analysis Method. We now demonstrate the 
numerical application of our approach over a finite-time interval for the case 7 = 3, the 
case also investigated in f^. We consider the set X = {(a;i,a;2) G M'^ : < xi < 2.0 
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Fig. 6.1. Sample trajectories of the DS system for 7 = 10.0. 



and < a;2 < 1.0} and check if the system (6.10 1, with 7 — 3.0, satisfies the 
conditions in Definition |4.1| for a finite-time uniform two-timescale set. Figure |6.2| 
shows the superposition of the forward and backward FTLEs, as functions of T, for 
a uniform grid of points in X; there is one fast-contracting exponent and one slow 
exponent uniformly on X. Taking /3 = 4.0, /i* = 1.0, /i-' = 3.0, n/f — 1, n'^ = 1, 



./« = 



0, To = 0.02, and Tc > 2.0, we find that the Definition 4.1 conditions are 



satisfied and we conclude that A" is a uniform two-timescale set. We note that with 
P = 4.0, Tc — 2.0 is sufficient. For the DS system, it can be verified that the 
timescale behavior is globally uniform, so that there is no upper limit on Tc- The 
FTLVs that approximate the fast and slow directions are if (T, x) and la^ (T, x) . The 
approximations _B-^'^(x) = span{lf{T,x)} and E''{x.) = sj)an{l^(T, x)} improve with 
increasing T sufficiently rapidly that, using T = 2.0, accurate approximations of the 
invariant subspaces i?^'^(x) and E''{x) can be achieved, and the invariant fast-slow 
splitting can be accurately approximated by 

TxK^ = E^'^iy^) ® ^'^(x) = span{l+(2.0,x)} ® span{l^(2.0,x)}. (6.15) 



Figure |6.3| shows points that are solutions to the orthogonality condition 
(f(x),lj^(T, x)) — for averaging time T — 2.0 illustrating that the FTLA method 
provides an accurate approximation of S. Both QR and SVD methods produced 
essentially identical results. 

In Fig. |6.4| the slow manifold tangent direction is compared to the approximation 
offered by the slow Lyapunov vector 1^ at points along the slow manifold S, indexed 
by xi. The Lyapunov vector 1^ provides a uniformly accurate approximation when 
a sufficiently large averaging time is used. T = 2.0 is large enough here, whereas 
T = 0.2 is not. 

6.2.2. Asymptotic Lyapunov Analysis. For the DS system, because the 
timescale structure is uniform on the entire state space (positive quadrant), the 
progress toward convergence in the first 2 units of time continues and it is possi- 
ble to compute the asymptotic Lyapunov exponents and vectors. The infinite-time 
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Fig. 6.2. Superposition of forward and backward FTLEs for the Davis-Skodje system illustrating 
uniformity. 



limits of the FTLEs can be determined analytically to be fj.^ = —7, and /ij = — 1. 
The backward time limits are (/u^,/i^) = (7, 1) = (— /i|J^, — /i;^). 

The fast and slow stable directions on the tangent plane at a point x = (xi, X2), 
are given by l^{T,x.) and 1^(T, x) respectively. We can analytically compute l2^(r, x) 
as the eigenvector of $(— T, x)-^<I>(— T, x) corresponding to the slow exponent in back- 
ward time, /i^(T, x). As time goes to negative infinity, l2^(T, x) can be shown to 
converge to 



I2 (x) = a{xi,X2) 



Xl + xi)' 



(6.16) 



where a(xi,X2) is a scalar function. For I2 to be a unit vector, a(xi,X2) should 



be chosen appropriately, 
converge to 



Similarly, as T goes to infinity, l^ (T, x) can be shown to 



l+(x) = 



(6.17) 



independent of x. 

If a point X is on S, then, using the asymptotic Lyapunov vector l^(x), the 
orthogonality condition characterizing points on S is in agreement with (6.11 1. These 



asymptotic results lend credence to the finite-time results, but the most important 
message is that in 2 units of time, the two-timescale behavior can be diagnosed and 
an accurate approximation of the slow manifold can be obtained. 
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Fig. 6.3. FTLA and ILDM approximations of the slow invariant manifold. 



6.2.3. 
Method). 



Invariant Slow Manifold Approximation Using Eigenvectors (ILDM 

The eigenvalues of Di in (6.141 are —7 and —1; in this case they indi- 
cate the two-timescale behavior correctly. Assuming that the eigenvector, denoted 
e*, associated with the slow eigenvalue —1, spans the slow subspace of the tangent 
plane, the ILDM method [IS] estimates points on S by computing solutions to the 
orthogonality condition < f(x), (e'*)-'- >= 0, where (e*)^ is orthogonal to e''. The 
slow eigenvector e'* can be obtained analytically and is 



e = 



lj_ (7+1) 



(7-1) 



Xi 



(6.18) 



The ILDM approximation to the slow manifold is 



X2 = 



Xi 



l + xi 



2^ 

1' 



(l-i)(l + xi)3 



(6.19) 



and is also shown in Fig. 6.3 The ILDM approximation is accurate around the equi- 
librium point (small xi) but gets worse away from the origin. The error is proportional 
to e^, where e = I/7, consistent with the analysis of the ILDM method in [21]. In 
Fig. |6.4| the slow manifold tangent direction, specifically the angle of the tangent rela- 
tive to the direction of the xi-axis, is compared to the approximations offered by the 
slow eigenvector e'* and the slow Lyapunov vector Ij^ at points along the slow manifold 
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Fig. 6.4. Comparing the slow manifold tangent direction to the directions of the slow eigenvector 
e" and the slow Lyapunov vector 1^ at points along the slow manifold S 



S, indexed by Xi. As expected, the slow eigenvector accurately approximates the slow 



manifold tangent direction near the equilibrium point (not shown in Fig. 6.4 1, but is 
in error farther away from the equilibrium point. 

6.3. 3D Nonlinear System: Hyperbolic Slow Manifold. Consider a non- 
linear time invariant (NTI) system 



Xi = 


XsXi 




X2 = 


= XfcX2 + ai{Xfc' 


- 2X,)xl 


±3 = 


= XfeXs + a2{Xfe 


- 2Xs)xl 



(6.20) 



For the numerical results, the constants are assigned the values A^ = —0.2, Xfc = 
—3, Aje = 3, and ai ~ a2 = 2. 

6.3.1. Finite-Time Lyapunov Analysis Method. First the FTLEs are com- 



6.5 



puted on a grid on the cubic region X = [—10, 10]'^ in the state space M.^. Figure 
shows a superposition of all the forward and backward FTLEs for different averaging 
times for the values of x on the X grid, indicating that X has uniform timescale 
structure. We see from Figure [675| that there is one fast stable exponent, one slow 
exponent, and one fast unstable exponent. Taking /3 = 6, ft" = 0.8, fi-^ = 3.5, 
n* = n^'^ = n^'^ — 1, To = 0.002, Tc = 3.0, the system satisfies the conditions given 
in Definition 4.1 to be a uniform two-timescale set. T = 3.0 will provide accurate 
FTLVs based on the bound given in Theorem 3.5. 

Having diagnosed two timescales and both fast contracting and fast expanding 
behavior, there may be a ID slow manifold and, if so, it is hyperbolic. Because 
there is sufficient averaging time for {E^{T,x.))^ — span{l]^(T, x),l;^(r, x)}, the ap- 



plication of the general result (5.11, to be a good approximation of the orthogonal 



complement to the corresponding invariant slow subspace, good approximations to 
the orthogonality conditions for points on slow manifold can be formulated. 

For the point x = [5 10 10]"^, Fig. 6.6 shows, as T increases, the distance, as 



defined in Proposition 3.2, between l^ (T, x) and Ij^ {T + AT, x) and l3(r, x) and 
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Fig. 6.5. Superposition of finite-time Lyapunov exponents in forward and backward time for 
grid on X. For forward time, the curves are red: fj,j{T,x), green: ^J{T,x), blue: fj,'^{T,x). For 
backward time, the curves are red: fi'^{T,x), green: fi2{T,x), blue: fi^(T,x), 



1^{T + AT, x), respectively, for AT = 0.006, confirming that the exponential bound 
given in Theorem 3.5 is satisfied. To compute the bound, we calculated the maximum 
eigenvalue of the symmetric part of the Jacobian matrix on the grid, and took Ti = 
0.006. The bound given in the theorem is conservative as stated earlier. In the plot, 
the bound has been reduced by a factor of four, to make it tighter. Although the 
FTLVs vary with x, the convergence behavior as T increases is represented uniformly 
over X by the behavior at [5 10 10]"^. The Lyapunov vectors for x = [5 10 10]"^ and 
T = 3.0 are 



If 



1^ 



1 



= [1.00 0.00 0.05]^ 
= [0.00 0.00 1.00]^ 
(6.21) 
Based on the results in (6.21 ) and the FTLVs at the other grid points of X, we chose 



= [0.00 - 1.00 0.00]^, 
= [1.00 0.05 0.00]^, 



= [0.05 0.00 - 1.00]^ 
= [0.05 - 1.00 0.00]^ 



Xi as the independent variable, because its coordinate axis is not parallel to any of 
the directions in {E^{T,x.))^. For each of the values on the Xi grid, we compute the 
values of X2 and x^ that satisfy the orthogonality conditions. The resulting finite-time 
approximation of the slow manifold for values of xi from -10 to 10 is plotted in Figure 
K7\ 

6.3.2. Exact Slow Manifold. For this problem, there is an independent means 
of determining the slow manifold, which allows the accuracy of the FTLA approach 
to be assessed. Over a time interval short relative to the fast timescale, yet long 
relative to the slow timescale, trajectories approach the 2D manifolds A4+ and A4^ , 
in forward (+) and backward (-) time respectively, given by 
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Fig. 6.6. Convergence of Lyapunov vectors, for x = [5 10 10]"' 




Fig. 6.7. Finite-time approximation of slow manifold 



M+ = {(xi, a;2, xs) € M^ | X2 + aixj = 0} 
M^ = {(xi, ^2, 2:3) e M^ I 2:3 + a2xl = 0} 



(6.22) 



The intersection of these sets is the slow manifold: that is S=A4^ f]Ai^ . These 
manifolds and their intersection are shown in Fig. |6.8| This intersecting of manifolds in 
the state-space (i.e., in the base space) is the counterpart of the subspace intersections 
in the tangent space employed in the FTLA method. 

At a point x G 5, the vectors normal to A4^ and A4^ are given by 771 = 
[2aiXi 1 0]^ and 772 == [2a2Xi 1]"^ respectively. Points on S, due to its invari- 
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(a) Manifolds: A^+(blue surface) and (b) Slow manifold: Intersection of 
M ~ (green surface) M ^ and M " 

Fig. 6.8. Slow manifold as intersection of manifolds 



ance with respect to the flow, satisfy the orthogonahty conditions 



= 



<r?i,f(x) > 

< [2aixi 1 0]^,f(x) > 



2aiXsxl + XfcX2 + ai(A/c - 


' 2X,)xl 


Xfc{x2 + aixj) 




<??2,f(x)> 




<[2a2Xi l]^,f(x) > 




2a2Xsxl + XfeX3 + a2{Xfe - 


- 2X,)xl 


Xfeix3 + a2xl) 





(6.23) 



where f (x) is the vector field given in ( 6.20 1 . We now show that the FTLVs conform to 
this geometry. As T ^ oo, l^(r, x) and 1 J (T, x) converge to tangent vectors of TxM~ 
and 1^(T, x) and l^{T,x) converge to tangent vectors of T^^A^"*". It follows that as T 
increases, 1^ should approach 772 and 1]~ should approach 771, as the results in (6.211 
indicate. For a given Xi, if we let {xi,X2,X3) denote the approximation of the slow 
manifold point {xi, 2x1, 2x1), then approximation error is [{x2+2xl)^ + {x3+2xl)^]^/^. 
We see that the error for FTLA approximation with T = 3.0, plotted in Fig. 6.9 



significantly smaller than that for the ILDM approximation. Figure 6.10 shows that 
the FTLVs closely approximate the normal vectors to the exact slow manifold for 
T>2. 

6.3.3. Invariant Slow Manifold Approximation Using Eigenvectors (ILDM 
Method). Applying the ILDM method, it is assumed that the eigenvector corre- 
sponding to the slow eigenvalue spans the slow subspace. The matrix for the linear 



variational equations corresponding to the system (6.201 is given by 

A, 

2ai{Xfc-2Xs)xi Xfc 
_2a2(A/e - 2Xs)xi A/e_ 



Dt 



the eigenvector corresponding to the slow eigenvalue. As, can be written as 

, . iT 

2^-)xi 



1 



-2ai(l 



'Af 



- 2^2(1 



%)^1 



Two linearly independent vectors orthogonal to Vj are 



(vi 



2aixi -4y^aia;i 



1 



(vi 



2aiXi 



Af^^aixi 



1 



(6.24) 



(6.25) 



(6.26) 
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Fig. 6.9. ILDM and FTLA slow manifold approximation errors 



Points on the slow manifold arc characterized as solutions to the orthogonality con- 
ditions 



((v^)i,f(x))=0, (K^)2,f(x))=0 



(6.27) 



For a given xi, the magnitudes of the errors in X2 and X3 relative to the correct 
values for S are 4ai(As/A/c)^a;f and 4a2(As/A/e)^a;f respectively. The slow manifold 
approximation error, as defined in section 6.3.3, for the ILDM method is plotted in 



Fig. 6.9 The ILDM error is similar to that for the FTLA method when the averaging 
time is short, but the FTLA method gives greater accuracy as the averaging time is 
increased. 

6.3.4. Roussel-Praser Method. The Roussel-Fraser method is based on pa- 
rameterizing the slow manifold as a graph and solving PDEs to obtain the graph. 
Parameterizing the slow manifold requires separating the state coordinates appropri- 
ately into independent and dependent variables for the graph. The FTLE/Vs provide 
a general means; however in this problem, the eigenvalues and eigenvectors suggest 
the appropriate choices. Assuming we have determined that xi is the independent 
variable and X2 and 2:3 are the dependent variables, we get the following PDEs 



ei 



£2 



dx2 
dxi 
dx3 
dxi 



X2 
Xi 

X3 

Xi 



+ {l-2ei)axl 
+ (1 - 2e2)axl 



(6.28) 
(6.29) 



where ei := X^/Xfc, and 62 :— Xs/Xfe- In general, the solution of the PDEs that arise 
in this method must be approached with numerical methods, and may be problematic 



LYAPUNOV EXPONENTS AND VECTORS FOR TIMESCALE ANALYSIS 



33 



0.8r 




dist( span{ri2(x)}, span{l3(T,x)} ) 



1.5 

Time 



2.5 



0.4r 

0.3 
0.2 
0.1 



dist( span{Ti^(x)},span{l^(T,x)} ) 



-2.5 



-2 




Fig. 6.10. Convergence of the FTLVs to directions normal to the slow invariant manifold 



even then [5 , but in this case the solution is straightforward. Assuming a;2(2;i 
Si^o^i^i ^^^ a;3(2;i) = J2'iLo^2^i^ ^^ obtain X2 " ' - — 



the correct graph for the slow invariant manifold. 



axl and 



-ax?. This is 



6.4. 5D System: Chaotic Attractor. The dynamics of a five-mode truncation 
of the Navier-Stokes equations for a two-dimensional incompressible fluid on a torus 
were studied by Frenceshini and Tebaldi [JT and later by Adrover et al. [T] . The five- 
dimensional dynamics depending on the constant parameter r, the Reynolds number, 
are 



Xi = 


-- -2a;i + 4a;2a;3 


X2 = 


: -9a;2 + 3a;ia;3, 


X3 = 


: -5x3 - 7a;ia;2 


X4 = 


= -5x4-XiX5, 


X5 = 


: — .T5 — 3X1X4. 



'4X4X5, 



(6.30) 



For r = 33, the case investigated here and in [1 , there is a chaotic attractor. 

We demonstrate the diagnostics derivable from FTLA, as a particular trajectory 
is followed in the state space. The trajectory begins at x = (1,10,1,10,10). The 
exponents corresponding to this trajectory, for averaging times T up to 15, are shown 
in Fig. |6.11| and were computed using the QR method. The FTLEs for T = 2.0 
were determined to be /i^ — —11.61, /i^ — —7.04, fi^ ~ —2.87, /i^ = —1.66, and 
HJ = 1.18; and the FTLEs for T = 15.0 were determined to be /i+ = -12.63, n^ = 
—8.27, n^ = —1.35, /i^ — —0.17, and /i^ — 0.42. For comparison, Adrover et al. pj, 
using an exterior-algebra based method '2<, determined the Lyapunov exponents to 
be nf = -11.97, ^J = -9.1, fi^ = -1.27, fi^ = 0, and ^i^ = 0.34, the sum of 
which is -22.0, for motion on the chaotic attractor. Their exponents are intended 
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to approximate the asymptotic Lyapunov exponents; for example, using their other 
computational method, they [T] use an averaging time of T=60. Because of our short 
averaging times, there should be no expectation that the FTLEs should match the 
asymptotic values. Also properties associated with the asymptotic exponents such as 
independence of the initial state x and a zero exponent associated with the direction of 
the vector field f(x) (if the attractor is a compact, invariant set without fixed points, 
e.g., a periodic orbit) will not in general apply to the FTLEs. 

In the basis provided by the backward Lyapunov vectors, we can write 

f(x) = E-=i"'.(x)i-(r,0(r,x)) (6.31) 

where wi, W2, W3, W4, and w^ are the coordinates of f (x) in the backward FTLV basis. 
The Euclidean distances from f (x) to the backward Lyapunov subspaces defined in 



(3.2) are 



rfi = dz,st(f(0(T,x)),/:2"(T,0(T,x))) = \wi\, 

rf2 = d*.st(f(0(T,x)),/:3-(T,0(T,x))) = {wl + wD^/^ 

d:,^dism<P{T,^)),Cl{T,<j>{T,^))) = {wl + wl + wlY'\ 



(6.32) 



We consider the system behavior along this trajectory over 2 units of time. The 
results in Fig. |6. 12] show how these distances evolve with time. The distances rfi, ^2, 
da, and d^, are analogous, but not equivalent, to the "mode amplitudes" in the CSP 
method [251 US) and in the work of Adrover et al. 11; we will refer to them as mode 



amplitudes anyway. The interpretation of Fig. 6.12 is not straightforward because (i) 
as time increases we are looking at what is going on in the tangent space at different 
points along the trajectory and (ii) as we move along the trajectory, the subspaces 
are converging to fixed subspaces due to the benefit of having averaged over a longer 
time interval. The mode amplitudes di, ^2, and ^3 of f (x) decay by 5 time constants 



respectively at the times 0.43, 0.71, and 3.00. The results shown in Fig. 6.12 indicate 
that the two fastest contracting "modes" have decayed significantly. However, for 
the decay of a particular mode to be relative to a fixed subspace, the subspace the 
mode amplitude is measured relative to must have converged (in a numerical sense, 5 
time constants being the criterion we are using). The convergence bound in Theorem 



3.5 



allows us to estimate that £^(T, ^(T,x)), C^{T,(l}{T,:>c)), and £4 (T, <?!)(r,x)) will 
converge in times of 1.09, 1.20, and 4.13, respectively, using 5 time constants for the 
appropriate A^ in each case. £2^(r, (/)(r, x)) and £^{T,(j){T,ii.)) converge sufficiently 
fast, as do the associated mode amplitudes, di and d2, that the decay of those modes 
is apparent in Fig. |6.12[ and there is a clear indication that the trajectory has reached 
a 3D manifold by T = 2 with tangent space given by C^{2,(f>{2,x)) at the base 
point 0(2, x). On the other hand, for 1^3, though the decay rate is large enough to 
be observed within 2 units of time, the subspace £4 (T, 0(T, x)) that it is measured 
relative to takes 4.13 units of time to converge, so it is not resolved by T = 2. 

Due to their convergence rates, for the second half of the time interval in Fig. |6.12| 
£2^(r, (/)(T, x)) and £3 (T, 0(T, x)) are very close to the fixed subspaces they would 
converge to if T were increased even farther and the timescale behavior were uniform. 
The FTLVs used to represent these subspaces were calculated by forward integration 
using the QR method with the initial state always x = (1, 10, 1, 10, 10). If we wanted 
better accuracy for the subspaces at trajectory points corresponding to the times 
between and 1, we would back up the initial condition along the trajectory under 
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consideration. For example to double T, one would start at 0(— T, x) and obtain 
£^(2T, 0(T, x)) by computing Af^{2T,(f>{~'T,x.)) with forward integration. 

This analysis for a single trajectory has allowed us to demonstrate the use of 
FTLA to produce the type of results sought by CSP and the method of Adrover et 
al. , and to illustrate the utility of the additional subspace convergence information. 
The computation of the FTLE/Vs for this example was more challenging. Whereas 
both the SVD and QR methods worked well and gave essentially identical results 
in the previous two examples, for this one the QR method permitted sufficiently 
large averaging time, whereas the SVD method did not. Further analysis of the 
system could involve computing FTLE/Vs more extensively and using the information 
to determine points on the lower dimensional manifolds to which trajectories are 
attracted, including the chaotic attractor. 
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Fig. 6.11. Finite-time Lyapunov exponents for the FT system with r = 33.0, for the initial 
condition x = (1, 10, 1, 10, 10). 



7. Conclusions. Two-timescale behavior of a finite dimensional, nonlinear time- 
invariant dynamical system on a not necessarily invariant subset of the state-space 
has been defined in terms of finite-time Lyapunov exponents and vectors, in a manner 
guided by the asymptotic theory of partially hyperbolic sets. Two-timescale behavior 
is characterized by a gap in the spectrum of finite-time Lyapunov exponents. There 
is a corresponding splitting of the tangent bundle into slow and fast subbundles de- 
fined by the finite-time Lyapunov vectors. The other desired property of a slow-fast 
splitting is that it is invariant under the linearized flow. In principle, determining an 
invariant slow-fast splitting requires computing asymptotic limits of the finite-time 
Lyapunov exponents. However, building on previous results, we have shown that 
under certain conditions the finite-time slow-fast splitting approaches an invariant 
slow-fast splitting exponentially fast as the time interval, over which the finite-time 
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Fig. 6.12. The distances between f{x) and the various finite-time Lyapunov subspaces. 



Lyapunov exponents and vectors are calculated, increases. The larger the spectral 
gap is, the faster the convergence. This is an important step toward establishing the 
feasibility of using finite-time Lyapunov exponents and vectors for timescale analysis. 
We have also provided evidence that the finite-time Lyapunov exponents and vectors 
more accurately characterize the timescales and associated geometric structure of the 
state-space than do the eigenvalues and eigenvectors associated with the "frozen-time" 
linear flow. 

When the tangent bundle has a slow-fast splitting, a slow manifold may exist. 
One approach for computing a slow manifold is to identify state-space points where 
the vector field is orthogonal to the directions normal to the slow subspace. In the 
intrinsic low-dimensional manifold method, the normal directions are calculated ap- 
proximately from the eigenvectors of the Jacobian matrix associated with the vector 
field, using the eigenvalues of this matrix to identify slow and fast directions. The 
alternative of determining the normal directions from finite-time Lyapunov exponents 
and vectors offers the potential for greater accuracy in determining a slow manifold. 
This advantage has been demonstrated in several application examples of increasing 
dimension and complexity. The examples illustrated that, consistent with existing 
theory, the accuracy of the eigenvector-based approach decreases as the curvature of 
the slow manifold increases and as the spectral gap decreases. The finite-time Lya- 
punov analysis method can yield more accurate normal directions even when there is 
significant curvature in the slow manifold. It can also yield accurate normal directions 
for a small spectral gap, if the gap is large enough relative to the available averaging 
time. 

Finite-time Lyapunov analysis of the linear variational equations provides an al- 
ternative diagnostic approach to eigen-analysis of the associated system matrix (the 
Jacobian matrix associated with the vector field). Though we have used this finite- 
time information to improve the performance of the intrinsic low-dimensional manifold 
type approach for determining points on a slow manifold, the finite-time information 
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could potentially be used (a) to suggest a transformation of coordinates leading to 
the standard form required for the analytical singular perturbation approach, (b) to 
initialize the basis vectors in the computational singular perturbation method, and 
(c) to guide the selection of independent and dependent variables in the application 
of the Roussel-Fraser partial differential equation approach. 

Further attention to numerical algorithms and additional application experience, 
along with direct comparison with other slow manifold determination methods, are 
needed to complete the development and assessment of the methodology. Although in 
this paper we have stopped with the computation of points on the slow manifold, our 
ultimate objective is to translate the geometric structure into reduced-order models. 
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